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ABSTRACT 

We examine the evolution of an almost circular Keplerian orbit interacting with unbound perturbers. 
We calculate the change in eccentricity and angular momentum that results from a single encounter, 
assuming the timescale for the interaction is shorter than the orbital period. The orbital perturbations 
are incorporated into a Boltzmann equation that allows for eccentricity dissipation. We present an 
analytic solution to the Boltzmann equation that describes the distribution of orbital eccentricity 
and relative inclination as a function of time. The eccentricity and inclination of the binary do not 
evolve according to a normal random walk but perform a Levy flight. The slope of the mass spectrum 
of perturbers dictates whether close gravitational scatterings are more important than distant tidal 
ones. When close scatterings are important, the mass spectrum sets the slope of the eccentricity and 
inclination distribution functions. We use this general framework to understand the eccentric i ties o f 
several Kuiper belt systems: Pluto, 2003 ELgi, and Eris. We use the model of Thol en et al.l ()2007l ) 
to separate the non-Keplerian components of the orbits of Pluto's outer moons Nix and Hydra from 
the motion excited by interactions with other Kuiper belt objects. Our distribution is consistent with 
the observations of Nix, Hydra, and the satellites of 2003 ELgiand Eris. We address applications of 
this work to objects outside of the solar system, such as extrasolar planets around their stars and 
millisecond pulsars. 

Subject headings: Kuiper Belt — planets and satellites — minor planets, asteroids 



L INTRODUCTION 

Several b i nary Kuiper belt objects (KBOs) have well-measured small orbital eccentricities (jNoll et al.l l2007f l. 
IStern et al.l (|2003f ) investigate numerically the forcing of the eccentricity of the Pluto-Charon orbit by interloping 
KBOs. They find that the system almost never possesses an eccentricity as high as the observed value of 0.003 
(jTholen et al.ll2007f l: depending on the model of tidal damping used, they find median values of 10"^ — 10~^. Our goal 
is to develop an analytic theory that describes the effects of a population of unbound perturbers on a binary orbit and 
can be applied simply to any binary, in the Kuiper belt or elsewhere. 

The interaction of a binary system with its environment has been studied extensively in the literature (jHeggid 
Il975t [Sigurdsson~fc Phinncv 1995; Yu 2002; Matsubayashi et al. 2007; Sesana et al. 2007). One interesting context is 
white dwarf-pulsar binaries, which are expected to be circular. For these objects pulse timing produces very accurate 
measurements of their orbi tal motion; s uch measuremen ts reveal that their eccentricities are typically very small but 
finite, around 10^"* — 10~^ ()Stairsll2dM ). iPhinneyl (|1992f ) investigated the effects of passing stars on the orbit of such a 
binary and found that for Galactic pulsars, the perturbations are sub-dominant compared to the effects of atmospheric 
fluctuations in the companion st ar. The higher density e nvir onment of a globular cl uster however can induce an order 
of magnitude higher eccentricity. iRasio fc Heggid (|1995D and iHeggie fc Rasiol (|1996D present a detailed account of the 
changes in orbital parameters for binaries in a stellar cluster. The work of these authors focuses on the regime where 
a perturbing body interacts with the binary on timescales longer than the orbital period of the binary. In the Kuiper 
belt, a single interaction between a binary and an unbound object occurs over a shorter timescale than the orbital 
period of the binary. We focus on this regime, where the perturbations to the orbital dynamics can be approximated 
as discrete impulses. 

The main result of this work is that we have identified the perturbative evolution of the eccentricity and relative 
inclination of a nearly circular binary orbit as a Levy flight, a specific type of random walk through phase space 
(jShlesinger et al.lll995| ). The entire distribution function of the eccentricity and inclination is then determined by 
calculating the frequency of perturbations as a function of their magnitude. We find a simple analytic expression for 
this distribution function. 

We take the following steps to arrive at our conclusion. In section [2] we calculate the effect of one perturber on 
a two-body orbit, examining separately the tidal effects of distant scatterings, close encounters with a single binary 
member, and direct collisions. We describe the effects of many such encounters in section [31 and write a Boltzmann 
equation that describes the distribution function of the orbital eccentricity and the inclination of the binary relative to 
its initial plane. The quantitative description of the binary's evolution given by this distribution function reveals its 
nature as a Levy flight. In section 31 we allow for a distribution of perturbing masses and discuss the different Levy 
distributions that result. 
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Fig. 1. — An illustration of the notation we use to denote the geometry of each perturbation. The dotted line is the almost circular orbit 
of the binary viewed at an angle. The dashed line is the path of the perturber, given by rp{t) = b + Vpt. 

We then use the analytic theory to examine the orbits of binary KBOs being perturbed by the other members of the 
Kuiper belt. Section [S] applies our analysis to several specific Kuiper belt binaries. We briefly discuss the relevance of 
this theory to other astrophysical systems in section [6l and summarize our conclusions in section [T] 



We use the following terminology to describe the geometry of the encounter between a single perturber and a two- 
body orbit. We refer to the two bound bodies as "the binary." The members of the binary have masses mi and m2, 
with a total mass labeled to;, = toi -I-TO2 and toi > TO2. The position of body 2 relative to body 1 is given by r^, and the 
relative velocity by v^. We distinguish between the magnitude and direction of a vector with the notation = rfcfb. 
We assume Vb ~ ^rti where D, is the orbital frequency of the binary. We write the orbital period as Torb — ^Tr/fl. 

We label the mass of the perturber nip. The position of the perturber as a function of time, Yp(t), is described by 
two vectors: rp{t) = b + Vpt. The vector b specifies the closest point of the perturber's trajectory to body 1, and Vp 
is velocity of the perturber relative to body 1. Each encounter geometry is uniquely specified by b and Vp under the 
constraint b • Vp = 0. Figure [1] depicts the arrangement of the vectors r^, v^, rp{t),h, and Vp. We assume Torb ^ b/vp 
so that we may ignore the motion of the binary during the interaction. We further assume that the effects of the 
gravity of the binary on the perturber are small; the perturber then travels along a straight path with a constant Vp. 
This assumption requires the criterion of 3> G{mb + mp)/b. If b is small, the perturber may collide with a member 
of the binary. In this case the assumption that the path of the perturber is unaffected by the gravity of the binary is 
true under the condition that Vp is much greater than the escape velocity of that member of the binary. The escape 
velocity from body 1 is defined v'^^^ ^ = 2Gmi/Ri, where i?i is the radius of body 1. 

We are assuming that the timescale of the interaction is much shorter than the orbital timcscalc, such that the 
perturbation instantaneously changes the velocities of the binary objects. The impulse provided to a specific member 
of the binary is found by integrating the acceleration caused by the perturber over its path: 



where the index j specifies whether the impulse Av^ and impact parameter hj are with respect to either the primary 
(j — 1) or the secondary (j — 2). For the primary, bi = b as we have defined it above. For encounters with the 
secondary, b2 is related to b by enforcing that it is also perpendicular to Vp. Thus we find b2 = b — + Vp{ri, ■ Vp). 

We consider the effects of such impulses on the full Laplace- Runge-Lenz vector, e = (v^ x JrVj/Gmi, — f, where 
H = Tb X Vb, the angular momentum per unit mass of the binary. The vector e has a magnitude equal to the 
eccentricity of the orbit, and points from body 1 towards the periapse. It responds to a small impulse Av according 
to the formula 



2. A SINGLE ENCOUNTER 




(1) 
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[2rf,(Av • Vb) - Vb(Av • rf,) - Av(rb • v^)] , 



(2) 



keeping terms up to linear order in Av. Since we have assumed the binary has very small eccentricity, the third term 
in equation [2] is negligible compared to the other two. 

The orbital plane of the binary is defined by the angular momentum vector H, and evolves according to AH = 
Ff, X Av. The impulses affect the direction of the angular momentum vector, and therefore alter the orientation of the 
orbital plane of the binary. We use the two-dimensional vector i to denote the components of H in the plane defined 
by the initial angular momentum. This vector, i, has a magnitude equal to sini, the sine of the inclination of the 
binary with respect to the initial orbital plane, and points from body 1 towards the longitude of the ascending node. 

The change in relative velocity given by a general gravitational scattering is given by Av = Av2 — Avi . The resulting 
change in the eccentricity vector is 



Ae = 2^^ 



In 



b2 ■ Vb _ b ■ Vb 
b2/rb b/n 



Vb 



b2 ■ h _ b ■ h 
b2/rb h/n 



The change in i is 



Ai = -2 



rup Vb 
nib Vp 



Vb 



b ■ n 



b2/rb b/n 



(3) 



(4) 



where n is the unit normal vector to the binary's orbital plane. For both the farthest perturbers and the closest, the 
dependence of equations [3] and [¥] on the impact parameter can be simplified. We discuss these limits in the following 
sections. 



2.1. Close Encounters 

Interactions with impact parameters greater than the radius of the primary or secondary but much less than the 
semi-major axis of the binary belong to what we call the "close-encounter regime." By definition the encounters in 
this regime of impact parameter are much closer to one member of the binary than the other. As a result the relative 
impulse experienced is dominated by the single impulse delivered to that body, |Av| « |Avj|. The changes in e and i 
are then given not by the difference of the impulses on each body, as in equations [3] and HI but by the effects of only 
the largest impulse. For the change in eccentricity we find. 



Ae = 2^^^ 
nib Vp b 



2fb{bj ■ Vb) - Vb{bj ■ h) 



and for the inclination. 



Ai 



^ mp Vb n 
nib Vp b 



Vb{bj ■ n) 



(5) 



(6) 



2.2. Distant Encounters 

For interactions where b ^ r;,, the impulse delivered to each member of the binary is almost the same. In this limit 
only the tidal difference in impulse affects the eccentricity of the binary. The perturbation delivered to the lowest 
order in rb/b is 



Ae - 2^J|^ (^) [h (4(r, • b){vb ■ b) + 2{fb ■ Vp){vb ■ Vp)) + h {l - 2{h ■ bf ~ (h ■ Vpf)] . (7) 

iPhinnevI {l992) der ives the special case of this formula for interactions that take place entirely in the plane of the 
binary. This formula is also equivalent to equation A24 of Heggic fc Rasio. (,1996i ). 
The change in i due to distant encounters is given by: 



mb Vp \ b 



Vb 



4(ff, • b){b ■ h) + 2{fb ■ Vp){vp ■ n) 



(8) 



2.3. Collisions 

Physical collisions between perturbers and body 1 or 2 cause the orbit to evolve impulsively. We define collisions to be 
any encounters where the impact parameter is smaller than the radius of the primary or secondary: 6 < ri or 62 < r2 . 
In this case the impulse is given by the conservation of linear momentum of the encounter: Av — x{rnp/mj)'Vp, where 
mj is the mass of the binary member involved in the collision (j = 1 or 2). The coefficient x accounts for the final 
momentum of the perturber. For an inelastic collision with nip ^ nij, x = 1- If the perturber is perfectly reflected, 
X = 2. The momentum loss from an impact crater can enhance this factor above 2 depending on the properties of 
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the colliding bodies (|Melosh et al.|[T99^ . For simplicity we assume that the mass of each binary member remains 
unchanged after each collision. 

The collisional impulse changes the eccentricity according to equation[2]and the orbital plane according to the change 
in angular momentum x Av. 

3. BOLTZMANN EQUATION 

The evolution of the eccentricity and inclination (relative to the initial orbital plane) is given by the sum of the 
perturbations the binary receives as it travels through a swarm of perturbers. From the average properties of the 
perturbing population, we can calculate a distribution function that describes the evolution of the orbit in a statistical 
sense. 

3.1. Eccentricity 

Since the perturbation in eccentricity is a two-dimensional vector, each component is added to the components of the 
existing eccentricity vector separately. As the binary experiences many perturbations, its eccentricity vector travels 
throughout this two-dimensional space. We write a distribution function /(e,i) that describes the probability that 
the binary will have an eccentricity in a small region (Pe. Assuming isotropic perturbations, there is no preferred 
longitude of periapse for the binary. It follows that /(e,t) = f{e,t) and the likelihood of finding the eccentricity in a 
small range de around e is 27re/(e, t)(ie. 

We define TZ{e') to be the frequency at which the binary experiences perturbations of magnitudes between e' and 
e' -|- de' . The frequency of perturbations with magnitudes on the order of | Ae| = e' is given schematically by e'TZ{e') ~ 
nvb'^, where n is the number density of the perturbers, v is the speed at which the binary encounters those perturbers, 
and b is the distance at which the binary encounters perturbers that cause a perturbation of strength e'. We make 
this calculation precise with the following integral: 

n{e')^ J S{\Ae{vp,h,m.p)\-e')J^{vp,m,p)vpS{h-Vp)d^hd^Vpdmp, (9) 

where !F{Vp,mp) is the phase space density per unit mass of the perturbers. The integral oi J-{^p,mp) over d^Vpdnip 
is the number density of the perturbers. We assume this density is uniform in the spatial dimensions and isotropic 
in velocity. It is normalized such that the total mass density of perturbers is given by p = J mpJ-{'Vp,mp)d^-Vpdrap. 
The factor of Vp in the integrand of equation [9] represents the velocity at which the binary encounters perturbers. The 
second delta function in equation [5] converts the volume element d^h to an element of cross-sectional area. The first 
delta function, i5(|Ae(vp, b, mp)| — e'), restricts the integral to include only the combinations of b, Vp, and mp that 
cause a |Ae| = e'. 

The evolution of the distribution function as a result of these perturbations is given by a Boltzmann equation that 
links the rate of change of /(e, t) to the interaction frequency. We write this equation as: 



df{e,t) 
dt 



= y'Me')[/(|e' + e|)-/(e)]dV (10) 

The function p(e') describes the frequency per unit of eccentricity space (d^e') at which a binary with eccentricity e 
is perturbed to the value e-|-e'. Since there is no preferred direction for the encounters, this function is axisymmetric, 
p{e') = p(e'). It is related to 7?.(e') by integrating over the angular direction of the phase space, 7?.(e') — J p{e')e'dio ~ 
27re'p(e'). 

We first derive p{e') for a simple scenario: a population of perturbers each with mass rUp and velocity Vp. To clarify 
this derivation, we present a qualitative treatment. The eccentricity excited by such a perturber with an impact 
parameter of order 6 >• rf, is about e' ~ {'mp/mi,){vb/vp){ri,/b)^ (Section m . Since the frequency of encounters with 
impact parameters b is proportional to 6^, and the size of the perturbation e' cx b~^, the frequency at which the binary 
is perturbed by an amount of order e' is therefore a power law: e'^p(e') oc e'^^. This power law is valid from very low e', 
caused by the farthest possible impulsive encounter, to e' ~ {mp/mb){vb/vp), the rare encounters with b ^ rb. We take 
into account the very rare occurrence of a physical collision, which excite eccentricities of order e' ^ {vrip / mj){vp / Vb) , 
in section [131 

Evaluating equation [5] using Ae(vp,b,TOp) given by equation [7] provides the exact form of p{e') for this scenario. 
We find: 

p{e') = McpTorb^' (11) 
47r e 

where Torb is the orbital period of the binary, and (Cg) — 1.89 is the average value of the angular terms of equation [7] 
(see Appendix). We note that the frequency of perturbations depends not on mp, but only on the total mass density 
of perturbers. It is also independent of Vp, as the lowered effectiveness of the faster perturbations is directly canceled 
by their higher frequ ency. These properties are typical of distant encounters with binaries, as evident in earlier work 
on binary dynamics (jBahcall et al. 1985). 

We can generalize eguation llOl bv including a term to account for dissipation of the binary's eccentricity: df{e, t) /dt — 
— div(/(e, i)e). We restrict our attention to mechanisms that reduce e at a timescale that is independent of e, 
e = — e/r^. The tidal dissipation of eccentricity obeys this form and is our main motivation for including such terms. 
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Since p{e') is a power-law, we can look for self-similar solutions to the time-dependent integro-differential Boltzmann 
equation, equation [TOl The frequency of perturbations p(e') does not depend on any special eccentricity, so the 
distribution function should depend only on the time t. We separate the distribution function into three parts: 
the time-dependent normalization, F{t), the time-independent shape of the function, g{x), and the time-dependent 
eccentricity scale, edt). These quantities obey the relation f{e,t) — F{t)g{e/ec(t))- We choose the normalization of 
g{x) such that / g{x)(Px = 1. We further choose that f{e,t) be normalized to 1 for all times; this constrains the 
normalization function to be F{t) — l/ec(t)^. 

Substituting f{e,t) = ec{t)^'^g{e/ec{t)) into equation \10\ we find two equations. The first specifies the time- 
independent shape of the distribution as a function of the dimensionless parameter x = e/ec{t): 

The solution to this equation has been presente d in several earlier works, where we inves tigate the eccentricity distri- 
bution of the oligarchs in a protoplanetary disk ([Collins fc Saril[2006l : ICollins et al.ll2007t l: 

g{x) = ^il + x')-'/^. (13) 

This function is the two-dimensional Cauchy distribution. The median and mode of this distribution are Xmcd = VS 
and Xmode = I/a/Z- The mean of this distribution is formally divergent; assuming there is a maximum value of x, 
Xu > 1, then ccmean ~ 2.3 log^Q (0.74a;,J . 
The eccentricity scale edt) is set by an ordinary differential equation, 

edt) = -edt)/Td + (Ce)Gprorb/2 (14) 

We note that Td and the terms on the right hand side of equation [Ml do not need to be constant in time; evolution of 
the binary {Torh{t))^ the perturbing swarm {p{t)), or the damping mechanism {Td{t)) can be treated by including the 
time-dependence of these quantities. 

We offer a reminder that edt) is the characteristic value of the entire distribution of eccentricity that the binary may 
attain. The probability is highest that the binary will have an eccentricity near the mode of the distribution, which is 
smaller than edt) by a factor of 0.7. The distribution is somewhat wide, and the confidence levels around the median 
value are large. The 66 percent confidence interval of x is 0.67 — 5.8, and the 95 percent interval is 0.23 — 40.0. 

Equal ions [T3l and [T4l present a new picture of the stochastic evolution of the binary's eccentricity. Often the evolution 
of a random variable is characterized by Brownian motion, in which the distribution of the random variable is set by 
the long term accumulation of many small perturbations. The typical value of such a variable grows as the square-root 
of time (written (x^) cx i^^^), and the probability of finding the system very far away from the typical value is 
exponentially low. The eccentricity of the binary evolves differently. The probability of finding the binary with an 
eccentricity larger than ec{t) only diminishes as a power law (equation ll3p . Physically, this refiects the probability that 
the binary received a single large perturbation to that state. The characteristic eccentricity, ~ ec{t) corresponds to the 
size of the perturbation that occurs with a frequency of about 1/t. The linear growth of edt) demonstrated by equation 
[141 reveals that the eccentricity of the binary does not reflect the accumulation of many small pertu rbations, but the 
single largest perturbation occurring in its history. This kind of random walk is called a "Levy flight" ([Shlesinger et al.l 
[19951 

3.2. Inclination 

The same analysis applies to the changes in angular momentum of the binary. Since |Ai| |Ae|, it follows that 
p(i') pie')- The evolution of inclination differs only in the coefficients that depend on the geometrical configuration 
of the encounter. The calculation of the coefficients is described in the Appendix. The self-similar distribution shape 
is a function of the dimensionless variable i/idt), where idt) is the time-dependent characteristic inclination. The 
following equation describes the evolution of idt)- 

Ut) = -^dt)/rd,^ + (C,)Gprorb/2 (15) 

where we have used Td^i to distinguish the timescale at which the inclination of the binary is damped, and (Ci) — 0.75, 
the average of the angular terms in equation [H The inclination is always measured relative to the orbital plane at 
t = 0- The distribution given by equation 1131 then describes the probability of the binary being inclined hy i = x idt) 
relative to its original orbital plane. 

4. A SPECTRUM OF COLLIDING PERTURBERS 

For many physical applications we must consider a range of perturbing masses and velocities and the effects of 
collisions onto the binary. In the single mass case discussed in section \JA[ the interaction frequency p{e') is set by 
the likelihood that the binary encounters a perturber at the impact parameter that causes such a change of e'. For 
perturbers that have different masses, the chance of experiencing a perturbation of magnitude e' depends on the 
combined likelihood that the perturber has the required impact parameter and the required mass to excite such a 
change. 
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To extend our analysis we set up several pieces of notation. Assuming that the mass and velocity distributions are 
independent, we consider T(rnp, Vp) — Tv{vp)Tm(jnp) . We restrict our analysis to velocity distributions with a charac- 
teristic value, Wo, such as a Gaussian distribution. We consider systems with differential mass spectra characterized by 
a power law: J-mi'fnp) cx rn.p'^ , valid from a minimum mass mmin to a maximum mmax- These functions are consistent 
with cond itions in the Kuipe r belt, where a power law mass spectrum and roughly Gaussian velocity spectrum are 
observed (|Luu fc Jewittl[200^ . We define the differential mass spectrum by 

^m{mp) (^0(7- l)/mo)(TOo/mp)^, (16) 

where ng is the number density of bodies larger than mass mg. In the literature the differential size spectrum of Kuiper 
belt objects is characterized as a power law in radius with index q; this is related to our index hy j — {q + 2)/3. In 
this section we discuss the p{e') and p{i') that result from several values of 7. 

4.1. 7 < 2 

The total mass density of perturbers for 7 < 2 is dominated by the perturbers with the largest mass, mmax- While 
perturbations of size e' are excited by all of the perturbers, the most likely perturber to cause a perturbation of this 
strength is the largest mass perturber. The dynamics of the binary are then the same as described in section [3. II with 
rrip = TOmax- The power law of p(e') oc e'~^ , based on distant encounters, is valid up to the eccentricity excited by a 
perturber of mass mmax interacting at a & ~ ri,, or for e' <C [m-n-iax/ 'rrn,){vh / Vp) (equation [7]). It is necessary only to 
know the total mass density p of the perturbing swarm in order to calculate the excitation frequency in this scenario, 
given by equation [TTl 

4.2. 7 = 2 

The power law 7 = 2 describes a special mass distribution where the frequency of encountering the few large 
perturbers at large impact parameters is the same as encountering the more abundant smaller perturbers at smaller 
impact parameters. Thus each logarithmic interval in impact parameter contributes the same amount to the frequency 
of perturbations by e', p(e'). The upper limit of impact parameters that can contribute to excitations of a given e', 
however, is given by the maximum mass perturber. The total range of contributing impact parameters then diminishes 
as e' approaches the eccentricity caused by the largest perturber interacting with b rt, e'^^^ = (TOinax/'7if,)(wf,/wo). 
Mathematically this behavior is determined by the integral of equation [51 which yields an excitation frequency of: 

Gnomoro.b log (2.1(e:„,,/e0) (Ce) .^„. 

for e' <C e^ax- The equivalent formula for the inclination excitations is: 

, ^ GnoTOoTorb log ((e^nax/^O) ,^gN 

For the smallest e' and j', the entire range of perturbing masses contributes to the interaction frequency. This occurs 
for excitations of the order {rritnin/ Trib) (vf, / vq) , below which the perturbation frequency is given by eauation llll 

4.3. 2 < 7 < 3 

The mass density of the perturbers when 2 < 7 < 3 is dominated by perturbers of the smallest mass, ?Ti,„in. Distant 
encounters by perturbers with this mass produce very small perturbations; for very low e' then, p{e') (x e'~'^, given by 
the simple model of section [3. II The upper limit of e' caused by these perturbers interacting with impact parameters 
6 ~ rf, is e' ~ imrain/mb){vb/vp). 

Perturbers with mmin cause eccentricity changes larger than this via close encounters, but these encounters are less 
frequent than interactions with perturbers of a higher mass and an impact parameter of order rb- Perturbations with 
a strength {niynin /'mb){vb/vp) >• e' >• (m,„ax/mb)(i'b/i'p) are most often excited by perturbers with impact parameters 
of ~ Tb and masses m ~ e'{vp/vb)mb. In other words the frequency of perturbations is directly proportional to the 
slope and normalization of the mass spectrum. 

In this case, the functions p(e') and p{i') cannot be determined using the simplifications to equation |3] afforded by 
very small or very large impact parameters. In general, the perturbation frequency for a mass spectrum of 2 < 7 < 3 
follows the power law p(e') oc e'~^'*'+^'. As an example we present the perturbation frequency for 7 = 25/12. This 
corresponds to q — 4.25, the best fit to observations of the Kuiper belt size distribution presented bv iFraser et al.l 
(|2008i) . We calculate from equation [21 

g/di-/!^ \mb uo/ 

It is simple to understand the relationship between equations II II and 1191 with the following argument. A perturbation 
of size e' that occurs via an interaction at a distance requires a perturber of mass about m' ~ e' (vo/vb)mb. If 
we interpret the total density in equation 1111 as only the density in bodies around m', then p' ~ m'J>n,(m') ~ 
'T'o(mo/m')^~^, and we recover the scaling of equation 1191 
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The integral over b and the angular variables of equation [4] yield a different coefficient for the perturbations to 
inclination: 



(20) 



, GnomoTorh f mo Vb\^^^^ 
P^^ ) = 7T77T2 — 

We relegate to the appendix the details of the integrals that produce the coefficients of equations [12] and [20] 

4.4. C'ollisional Perturbations 

The integral of equation [9] over impact parameters from to rj produces the frequency of perturbations to the binary 
by collisions on member j. Since the size of the impulse from a collision does not depend on the impact parameter, 
it is the mass of the perturber that dictates the size of the eccentricity perturbation. Accordingly, the frequency of 
perturbations as a function of e' reflects the frequency of collisions as a function of rup. The frequency of coUisional 
perturbations does not depend on Wmax or mmin regardless of the slope. However, the limits of the mass distribution 
specify the lowest and highest perturbations achievable via collisions: x{'mmin/fnj)(vo/vb) < e' < xim-ma.^/ rn-j) (vq / Vb) . 
In this range of e', for any value of 7, the perturbation frequency due to collisions is 



Pie') = 



e/7+1 



mo 



7-1 



— I ( - I ( ^ I K 



(7-i)(jr^) 

2tt 



(21) 



where {DJ^^) is the average of the angular dependence of Ae from collisions to the power of 7 — 1, and Ky = 
Vq'^ J v'^^'^J-'y(vp)dvp. If J^v{vp) is proportional to a delta function, S{vp — vq), then V-y = 1 for all 7. If the velocity 
spectrum were Gaussian, such that Tv{vp) oc exp(— (wp/wo)^), then V-y — 2r((3 + 7)/2)/y^. The frequency of 
perturbations to the relative inclination by collisions is the same as equation I2H replacing the integrated coefhcient 
{D2~^) with the appropriate calculation made from the coefficients of |Ai|. 

Although we use rj to represent either member of the binary, it is clear from equation 1211 that the collisions onto 
the smallest body have the largest effect on the orbit. The ratio of the perturbation frequency through collisions, 
p(e')coUisions fcquation 12 ip to the frequency of gravitational scatterings, p(e')gi.avity (equation I19p . is, for mass distri- 
butions of 2 < 7 < 3, 



P(e ) collisions — Q f ^ J ^ 
P(e') gravity ' \rb J 



nib I 

X — — 

Vb 



(22) 



where we have evaluated the coefficients for 7 — 25/12. The choice of 7 does not change these coefficients dramatically. 

4.5. Eccentricity Distributions 

The distribution given by equations [T3] and [TJ] were derived in the context of p(e') cx e'~^. As long as p(e') 
follows a power law with e', we can write a self-similar distribution function f(e,t). We write a generic function, 
p{e') = Ppe' '•^^'''1 to account for the different slopes caused by different mass distributions (for 3 > 7 > 2, r/ = 7; 
for 7 < 2, ?7 = 2). The derivation of the distribution function proceeds analogously as in section [3.11 Equation 
[TUI becomes two equations: a dimensionless integro-differential equation that specifies the shape, and an ordinary 
differential equation to specify the evolution of the eccentricity scale edt). The general version of equation [Ml is: 



S)/Td + 27TPo/e,{t)'J- 



,{t) cx 



(23) 

For all of the p(e') 



In the limit of no eccentricity dissipation (r^ — > 00), equation [^ shows that Cc 
discussed in section [4] the growth of edt) is always faster than i^/^. 

The shape of the distribution function is determined through a Fourier transform of the gener al version of equation 
[H For slopes of 1< < 3, g{x) = J cos(k • x) exp(-|k|''-i)d2k (ISatolll999l: [Collins et al.ll200l . While there is only 
a closed form solution for rj — 2, given by equation [T51 all of these functions are flat at low x and fall off like x~^'^~^^\ 
In fact, it is easy to show from eq nation 1 1 01 that the high e tail is given by 



/(e»e,(t)) =p(e)t/(7-l), (24) 

when t <C t^. For equilibrium distributions where edt) = 0, t is replaced with r^, the timescale for the dissipation. 

When p{e') cx e'^* or steeper, the accumulation of the smallest perturbations over time is more effective at raising 
the eccentricity of the binary than single large perturbations. In this case, the evolution of the eccentricity follows 
standard Brownian motion, where the distribution function is a Gaussian, and edt) cx t^^^. 

5. KUIPER BELT BINARIES 

In this section we compute e^t) and ic{t) for several Kuiper belt binaries. The "binary" of section [5] now refers to 
a bound pair of Kuiper belt objects, and the "perturbers" are all of the other members of the Kuiper belt. 
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For the highest mass KBOs, the size spectrum is weh determined to be a power law with an index sHghtly greater 
than q — A. The lowest mass bodies, of about 30 km in radius, are less frequent than predicted by a single power law, 
however the pa rameters of a more genera l model are still under inv estigation (jTruiillo k, Brownll2001l : iLuu &: Jewittj 
I2002t iPan fc S ari 2005; iFraser et al]|2008t iFuentes fc HolmanI [20081 ) . For this section we use the best fit of a single 



power law model to the high mass part of the spectrum provided bv iFraser et alJ ((2008') , who find q = 4.25 and a 
number density of 1 body per square degree brighter than magnitude 23.4. We assume an average distance of 40 AU 
to the Kuiper belt and a depth of 20 AU to find a volumetric number density no = 3 x 10^*^ cm^'^. To convert 
the magnitudes of the objects to physical sizes, we assume a constant geometric albedo of 0.04, a constant physical 
density of 1 g cm^^, and take the R-band apparent magnitude of the Sun to be -27.6. We find that the magnitude 23.4 
corresponds to a mass mo = f .75 x 10^^ g, equivalent to a radius of 75 km. Most of the objects found between 30-50 
AU are inclined by about 5-15 degrees relative to the plane of the solar system, and have heliocentric eccentricities of 
0.1-0.2. 

5.1. Perturbations by a Disk 

Our analysis so far has treated the perturbing bodies as unbound objects moving relative to the binary with a 
constant velocity. When the perturbers are part of a disk orbiting the central star, the orbital elements of the disk set 
the parameters of the perturbation frequencies we calculate in section [3l 

The relative velocity between KBOs, when they interact, is set by the size of their eccentricities and inclinations, 
Vp end^H, where the subscript "H" denotes a heliocentric orbital quantity. We assume a constant perturbing 
velocity with Vp = 1 km/s, which corresponds to the typical heliocentric eccentricities and inclinations of KBOs. We 
assume that these encounters occur isotropically in the frame of a binary, however this is not accurate. A more detailed 
calculation of the angular distribution of relative velocities will only affect the coefficients of the perturbations. The disk 
does not specify a special direction for the perturbation vector Ae, so the perturbing frequency and the distribution 
function retain their axisymmetry. The influence of the central star on the binary and the perturbers adds another 
constraint to our assumption of impulsive encounters: the timescale for an interaction must be shorter than the orbital 
period around the star: b/vp <C 1/Ojj, or equivalently, b <^ ena. This guarantees that the relative velocity is constant 
during the interaction. 

If the orbit of the binary is much different than the typical KBO orbit, there are several modifications to perturbation 
frequencies experienced by the binary. One modification is due to the finite height of the disk of perturbers. This 
height is set by their inclinations around the central star; for the Kuiper belt we refer to the average inclination as 
{i)KB- A binary with heliocentric inclination icoM {'i)kb never travels above or below the perturbing disk height 
and therefore experiences the maximal frequency of perturbations. If icoM ^ (*)a's, the binary spends most of its 
orbit outside of the perturbing swarm. The frequency of perturbations to such a binary is reduced by the fraction of 
the time the binary leaves the disk, proportional to {i) kb I icoM- The eccentricity of the binary in the disk reduces the 
effective density of perturbers in a similar manner if the epicycle of the binary carries it outside of the region populated 
by perturbers. 

If the heliocentric eccentricity or inclination of the binary is much greater than the typical values for the Kuiper 
belt, the relative velocity between the binary and a perturber is primarily due to the non-circular heliocentric motion 
of the binary. Gravitational interactions depend weakly on vq so their frequency does not change much in this case. 
Perturbations by collisions, however, become more important if vq is increased due to this effect (equation [22]) • 

5.2. Pluto et. al. 

Pluto is the second largest known Kuiper belt object, with a radius of about 1100 km. It has a semi-major axis of 
39.5 AU and its orbit is inclined relative to the ecliptic by 17°. Its largest satellite, Charon, contains abo ut one tenth of 
the t otal mass of the system. Recent observations have revealed two smaller satellites. Nix and Hydra ([Weaver et al.l 
l2006l ) . These satellites have small eccentri cities and are roughly co-planar with Charon. Numerical simulations of 
collisions between similarly sized objects by ICanupl ([20051 ) produce binaries with orbits similar to Pluto and Charon. 
The circularity and co-planarity of Nix and Hydra lend additional weight to a coUisional origin of the system. 

The triple system of Pluto and its moons is a valuable test case for the dynamics we have presented. For an isolated 
binary it is impossible to know the initial orbital plane. The relative inclinations of the moons of Pluto can be measured 
directly assuming their formation was co-planar. Furthermore, the perturbing swarm for all three Pluto-moon pairs 
is the same. A major issue in comparing our analytic calculations to the observations is that the large mass ratio of 
Charon to Pluto causes significant non-Keplerian effects in the orbits of the outer satellites. We first re-examine the 
published observational model of their orbits to separate the relevant motion of the outer satellites from the forced 
motion due to Charon. We then compare the resulting eccentricity with our predicted values. 

5.2.1. Orbital Model of Tholen et al 

A model of the observations of the Pluto system has been presented bv lTholen et al.l ([2007f ). who fit the parameters 
of a four-body numerical integration such that the simulation agrees with the observations. Such work is necessary, 
as it has been show n that the observations cannot be consistently modeled by three non-interacting two-body orbits 
([Weaver et al.ll2006D. 



The model of iTholen et al.l (|2007( ) presents a full set of osculating el ements describ i ng th e orbits of Charon, Nix, and 



Hydra. The orbit of Charon is virtually unaffected by Nix and Hydra; iTholen et al.l ([2007[ ) measure the eccentricity of 
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Time (days) 

Fig. 2. — The distance of Nix (lower panel) and Hydra (up per panel) from the Pluto-Charon barycenter, in units of Pluto radii, as a 
function of time, in an integration of the parameters found bv lTholen et al] 120071) . Nix and Hydra are treated as massless test particles. 
The origin of the time coordinate is arbitrary. 

Charon to be 3.48 ± 0.04 x 10^'', and the period of its orbit is 6.387 days. Since the combined potential o f Charon and 
Pluto is significantly non-Keplerian, the elements of Nix and Hydra vary significantly during their orbits. iTholen et al.l 
()2007l) average the osculating semi-major axis to find an orbital period for these satellites of 25.49 days and 38.73 days 
for Nix and Hydra respectively. The osculating eccentricities of Nix and Hydra both oscillate between zero and about 
0.2; for each satellite oscillations at the frequencies of its own orbit and that of Charon are visible (their figure 4). The 
orbital planes of the satellites relative to Charon's are tilted by 0.15 degrees for Nix and 0.18 degrees for Hydra. Each 
plane precesses relative to the plane of Charon, however the angle of the offset remains constant. 

5.2.2. A Different Interpretation 

For two body motion, the Keplerian elements are constant and indicate the shape of the orbit in space. Osculating 
elements that describe motion in significantly non-Keplerian potentials, such as the combined potential of Pluto and 
Charon, may vary on timescales shorter than the orbital period of the satellite. When this is true, relating the 
osculating elements to t he shape of the orbit can be misleading. The average value of the osculating eccentricity of Nix 
is 0.015 in the model of iTholen et al.l (|2007f ). however the motion of Nix relative to Pluto never resembles an ellipse 
with such an eccentricity. 

We re-examine the model provided bv lTholen et al.l (|2007f ) by reproducing the numerical integration based on the 
Pluto-centric positions and velocities of Charon, Nix, and Hydra published in their table 1. We set the masses of 
Nix and Hydra to zero to eliminate their secula r inter actions with each other. Instead of examining the osculating 
elements, we adopt the approach of iLee fc Peald (|2006f ) and characterize the orbits of Nix and Hydra based on their 
position as a function of time from the Pluto-Charon barycenter, plotted in figure [2l The units of distance are Pluto 
radii, defined as Rp — 1147 km. 

Although short oscillations on the timescale of Charon are visible in the top panel of figure [H they are very small 
compared to the oscillations that occur on the timescale of Hydra's orbital period. To parametrize Hydra's orbit we fit 
the function ro{l + e cos{Kit + wi)) to the first 200 days of the numerical model. Because for a non-Keplerian potential 
the radial epicyclic frequency differs from the orbital frequency, we calculate the average angular frequency by fitting 
a straight line to the angular position of Hydra as a function of time, f{t) = flit + Aq. The results are written in table 
[TJ We interpret ei as the orbital degree of freedom in the combined potential of Pluto and Charon that is analogous 
to the eccentricity of a two-body orbit. 

The motion of Nix (bottom panel of figure[2]) appears more irregular than that of Hydra. We find the position of Nix to 
be well-described by a model of three epicycles with different frequencies: r{t) — ^0(1 + J2k=i 2 cos{Kkt + iL>k)). The 
best fit values are printed in table 1. We distinguish the cause of each epicycle by its period. The combined potential 
of Pluto and Charon oscillates with frequency of rJcharon ~ ^^Nix! motion being forced by this potential should occur on 
integer multiples of this frequency. Using the numbers in table[Tl we see that 27r/ (r2Nix + '«2) = 27r/ (ilNix + K3/2) — 6.39 
days. The second and third epicycles in our fit correspond to motion at the first and second harmonic of Nix's relative 
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TABLE 1 

Best fit values to the epicyclic models of the radial motion of Nix and Hydra. 



ro/Rp 


ei 


2-n 1 K\ 


£2 




£3 


27r/K3 


27r/ni 






(days) 


xlO-3 


(days) 


xlO-3 


(days) 


(days) 


Nix 46.805(5) 


2.96(3) 


25.22(2) 


1.25(3) 


8.599(8) 


1.38(3) 


4.298(1) 


24.8505(5) 


Hydra 62.237(1) 


5.595(2) 


38.535(15) 










38.20(1) 



Note. — The motion of Nix is fit with three epicyclic terms, while the motion of Hydra is only 
fit with one. The parenthesis indicate the 95 % confidence level of the fit around the last digits. 

orbital frequency. We therefore interpret the first term, with a size of ei = 3 x 10"'^ and a period close to Nix's orbital 
period, as analogous to the two-body eccentricity. 

We perform another integration of the best fit initial conditio ns from|Tholen et al.l ()2007f ) to investigate the secular 
effects between Nix and Hydra. We use the best fit masses from lTholen et al] l 2007f ) for the two outer satellites. Since 



the motion of Hydra is dominated by a single epicyclic frequency, the variation in the size of its epicycle is apparent on 
the timescale of several years. To determine the effect of secular variations on Nix, we fit the same three-component 
epicyclic model to five orbits at t ~ 5 years. In the best-fit model to these later orbits, the only difference compared 
to the model of table [T] is in ei, the epicycle with a frequency close to Hydra's orbital frequency. This is further 
confirmation that the degree of freedom represented by ei is analogous to the two-body eccentricity. 

5.2.3. Theoretical Distribution 

To compute the distribution of eccentricities and inclinations expected of Pluto's moons, we solve equation l23l for each 
of the moons, given the interaction frequencies specified by equations [TOl and The only remaining parameters to 
evaluate are the damping timescales for the eccentricity and inclinations of each satellite. We use the standard formula 
for the d amping of eccentricity due to the tidal force o f the primary acting on a secondary that is in synchronous 
rotation (|Yoder fc Peald[l98ll : lM7irrav & Dermottlll999f ): 

r.. = i.Q.(l + A.)^(^)'^, (25) 
00 mi \r2 J \l 

where Q2 is the dissipation function of the secondary, and jl2 = 19/ir2/ {2pGm2) is its effective rigidity, a ratio between 
the material strength of the secondary and its self-gravity. The damping rate of eccentricity due to tides of the primary 
acting on the secondary, Td,i, if the primary is also rotating synchronously with the orbit of the satellite, is given by 
equation 1251 with the quantities specific to the primary switched with those of the secondary and vice versa. 

Pluto and Charon are known to be in a double-synchronous state of rotation, where the spin period of each body 
is equal to the 6.4 day orbital period. In many binaries, only the spin of the secondary is synchronous with the 
orbital period. Tides on the primary then raise the eccentricity. Double-synchronous systems, however, experience 
damping due to both the tides on the secondary and those on the primary. Assuming a water-ice composition for Pluto 
(/X = 4 X 10^*^ dynes cm~^), we calculate the eccentricity damping timescale due to tides raised by Charon, Td,i from 
equation[25lto be 5.1 Myrs. The shortest damping timescale due to tides from Pluto acting on Charon, Td.2 is found by 
assuming Charon is also made of water-ice; we find in this case a timescale of 8.2 Myr. The longest timescale assumes 
a rocky composition (/i = 6.5 x 10^^ dynes cm~^); we find this corresponds to 133 Myr. The overall damping of the 
system is given by the sum of the damping rates. The short damping timescale of tides on Pluto prevents Charon 
from contributing significantly to the combined effect of both tides, reducing the importance of its composition. The 
longest eccentricity damping timescale that results from both tides is 4.9 Myr. The inclinations of the outer satellites 
relative to the Pluto-Charon plane are also damped by tidal dissipation. For a circular synchronous orbit the timescale 
for inclination damping is longer than the timescale for eccentricity damping by a factor of ^ i^^. We ignore the 
damping of inclin ations in equa tion 1151 for all three satellites. 

As discussed in iTholen et al.l (|2007 ) and section 15.2.21 secular interactions between the satellites are visible in the 
long term calculations of their orbits. For the best-fit values of the masses of Nix and Hydra, their eccentricities 
are modulated on the order of 10% over timescales of years; we neglect these fluctuations for this work. It is more 
important in this model to determine whether secular evolution can cause the eccentricity of Nix or Hydra dissipate 
via Charon's orbit. 

We use line ar secular theory to desc ribe the coupled evolution of the eccentricity and longitude of periapse of 
each satellite (jMurrav fc DermotHllQQQf ). We find that the undamped secular evolution agrees qualitatively with the 
numerical orbit determinations. We add a term to the differential equations describing Charon's eccentricity that 
reduces it at a constant timescale (echaron = — echaron/Td)- The frequencies of the oscillations of the eigenmodes of the 
solution are practically unchanged by this term, however each eigenmode gains a dissipative factor. Quantitatively, 
only one eigenmode is damped on timescales shorter than than 4.5 Gyr. By integrating the damped secular equations 
with different initial periapses, we determined that the secular interactions do not cause substantial damping of Nix 
and Hydra. 

Equation [22] gives the frequency of perturbations due to collisions of perturbers onto each moon relative to the 
frequency of perturbations caused by gravitational scattering, equation 1191 For Charon, the coUisional perturbations 
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increase p{e) by only 2 percent. Since Nix and Hydra are smaller, perturbations by collisions have a greater relative 
effect; however it is only a 20 percent contribution to the total perturbation frequency for Nix and 15 percent for 
Hydra. We solve equation!^ to find edt) and ic{t) for each of Pluto's moons. 

For Charon we find Cc = 2.6 x 10~^, and ic = 0.029°. This value of Cc corresponds to the mos t likely perturbatio n 
during a damping timescale of 4.9 Myr, and is much smaller than the observed value of 3.5 x 10"'^ (jTholen et al.l[2007b . 
Using equation [24l we calculate that given this value of Cc, the probability of Charon's eccentricity being as high as 
its observed value is 0.2 percent. 

For Nix we calculate ec(4.5 Gyr) = 4.8 x 10"^ andjc(4.5 Gyr) = 0.1°, and for Hydra, 7.1 x 10"^ and 0.15° respectively. 
The distributions specified by these values are quite consistent with the free eccentricity we determine in table [TJ 

5.3. Other Interesting KBOs 

Two other Kuiper belt objects have satellites on low eccentricity orbits: 2003 ELgi, and Eris. Along with Pluto these 
are three of the four most massive KBOs known, all with radii of about 1000 km. 2003 ELsihas two k nown satellites. 
The largest has a 50 day orbit and a measured orbital eccentricity of 0.050 ±0.003 ([Brown et al.ll2005h . An additional 
smaller satellite orbits 2003 ELgi with a period of about 35 days (|Brown et al.ll2006h . The orbital parameters of the 
inner satellite are unconstrained, however the relative inclination between the two is about 40°. The masses of the 
s a.tellites are n e gligib le compared to the mass of 2003 ELgi. The heliocentric inclination of the system is 28°. 

iBrown et all (|2005f) argue that if the tidal response of 2003 ELgi and its large satellite are fluid-like, tidal interactions 
should damp their eccentricity on a timescale of about 300 Myr. With these parameters we use equation[23lto calculate 
an equilibrium Cc — 4.3 x 10"**. The distribution with this eccentricity scale predicts an observed eccentricity of 0.05 
at a probability of three percent. However, for smaller bodies, internal elastic forces dominate the tidal deformation 
of their shape; it is more reasonable to assume that the tidal response of the satellite is characterized by its material 
strength. Then, the tides raised on the primary have the greatest effect and the eccentricity of the system grows on 
the same timescale as the growth of the semi-major axis. Forced eccentricity growth and an evolving orbital period 
can be incorporated into equation 1231 However, these corrections are only an order unity correction since the growth 
timescale, by definition, is comparable to the age of the system. Assuming Torb is fixed and ignoring the eccentricity 
growth, we calculate ec(4.5 Gyr) — 0.0052. The 95 percent confidence interval around this ec is 0.001-0.2; the observed 
eccentricity of 2003 ELgiis within this range. 

The dwarf planet Eris is orbited by the satellite Dysnomia. Observations have shown an upper limit to their 
eccentricity of 0.013 f Brown et al.ll2006l ). The system has a 15 day orbital period, and orbits the sun at a semi-major 
axis of 67.7 AU with an eccentricity of 0.44 and a heliocentric inclination of 44°. In addition to the reduction in 
effective perturbing density caused by the inclination, the high eccentricity reduces the effective perturber density by 
an additional factor of 0.09. The semi-major axis of the binary is consistent with 4.5 Gyr of tidal evolution away from 
an initially very close orbit; if the tidal response of the secondary that of a strength-less fluid, then its eccentricity is 
damped on a timescale of 50 Myr. These parameters yield an Cc — 2.2 x 10~^. However, if the material strength of 
the secondary is stronger than its own self-gravity, then the tides raised on the primary cause the eccentricity of the 
satellite to grow. In this case the relevant timescale is the age of the system, and we find that ec(4.5 Gyr) = 1.0 x 10^*. 
Both values are below the observed upper limit. 

In addition to the high mass ratio and low eccentricity Kuiper belt binaries, there are other known binaries of almost 
equal mass on moderately eccentric orbits. The binary 1998 WW31 is an example of such an ob ject: both me mbers 
have a radius of about 50 km, an orbital period of 574 days, and a mutual eccentricity is 0.817 (|Veillet et al.l [2002^. 
Even though our analysis is derived in the low eccentricity limit, we can use equation [23] to estimate approximately the 
eccentricity expected from impulsive encounters; we find ec(4.5 Gyr) = 0.31. This moderate characteristic eccentricity 
is consistent with the high observed value. Other binaries with orbital periods on the order of a year will have acquired 
large eccentricities through their interactions with the other Kuiper belt objects. 

6. OTHER BINARY SYSTEMS 

Our analysis holds for any two-body orbit perturbed isotropically in the impulsive limit. As binary orbits are 
prevalent in astrophysics, we briefly discuss several other examples. 

The asteroid belt harbors many binaries with well determined eccentricities. The mass spectrum of the asteroid belt, 
however, is much shallower than that of the Kuiper belt: the largest asteroid, Ceres, contains a third of the total mass 
of all asteroids. A binary asteroid is then perturbed mostly by the largest objects that it encounters. To calculate 
p(e') accurately, it is necessary to model the neighborhood of that binary. The asteroid belt is also coUisionally active 
so its binaries may not be coeval with the whole solar system. We postpone a detailed analysis of the binary asteroid 
population for a future work. 

A well-measured class of binaries outside the solar system are millisecond pulsars with white dwarf companions. The 
tidal damping between the pulsar and its companion in the phase before the companion becomes a white dwarf is very 
short, indicating tha t during this phase the eccentricity of the binary should b e smaller than the observed values of 
around 10~^ — 10^^ ([Stairsll2d0^ . To explain the observations. iPhinnevI ([19921 ) presents the following model. As the 
companion star becomes a white dwarf, random fluctuations in the atmosphere of the star cause irregular motion in 
the orbit of the binary. These motions are reflected by a small eccentricity that remains s ince the t idal in teractions 
between the white dwarf and the neutron star cannot damp the system. The model of iPhinnevI ([19921 ) produces 
eccentricities for these systems that match the observations well. 

These binaries are perturbed by encounters with other stars in the galaxy; we can calculate the contribution to 
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their eccentricities by the distant stellar interactions. The perturbation of these systems by other stars falls into the 
simple regim e of only distant interacti ons described in section [3TT1 A typical volumetric mass density for field stars is 
O.IMq pc~'^ (|Holmberg fc FlY"iinll200ClD . Given this density, we calculate the characteristic eccentricity of these systems 
to be 

Typical orbital periods are between 1 and 10 days, and the ages of these systems are on the o r der of Gyrs. We find 
then that edt) is several orders of magnitude lower than the observed eccentricities. iPhinnevI ()1992[ ) also concludes 
that the perturbations from other stars cannot be responsible for the eccentricities of the binary pulsars. Since we have 
calculated the distribution, however, we can estimate more accurately the likelihood of achieving these eccentricities 
by only distant stellar perturbations: less than 0.1 percent. 

Globular clusters can have densities many orders of magnitudes higher than the average galactic density, such that 
distant perturbations to the binaries may be important. However, in a cluster the interactions between a binary 
and a star are not typically in the impulsive interaction regime. Instead the orbits of the perturbers are affected by 
the gravity of the binary, and the interactions occur over several orbital periods. Analytic work on the eccentricity 
perturbations in this regime has been performed by Rasio & Hcggic (1995) and Heggie & Rasio (1996). 

The characteristic eccentricity caused by distant stellar passages on the orbits of extra-solar planets is also given by 
eauation l26l These eccentricities are too low to be reflected in the current sample of known extra-solar planets. As with 
the pulsar binaries, the distant interactions may play a role in setting the eccentricity distributi on of long period planets 
found in a dense stellar cluster. For most extra-solar planets however, planet-disk interactions (jGoldreich fc Sarill2003| ) 
or planet-planet scatterings (,Rasio fc Ford, 199 6) are probably the source of their eccentricity. 

7. CONCLUSIONS 

We have calculated the effects of impulsive perturbations and collisions on a nearly circular Keplerian orbit. If the 
swarm of perturbers encounter the binary isotropically in space, we can write a distribution function that describes 
the probability density for the binary to have a given eccentricity or inclination relative to its initial plane. The 
growth rate of the binary's likeliest eccentricity and inclination depends on the mass spectrum of the perturbers. For 
shallow mass distributions {q < 4) it is the distant encounters that set the binary's eccentricity and only the total 
mass density of perturbers is important to the evolution of the binary. For steeper mass distributions of g = 4 — 7, it 
is the interactions at about the semi-major axis of the binary that dominate the frequency of perturbations. Only the 
normalization and slope of the mass spectrum set the distribution of eccentricities in this regime. 

The assumptions of this model are valid in the Kuiper belt. Our calculations match the observations of Nix and 
Hydra very well. For Eris and 2003 ELgi, the observations lie within the 95 percent confidence intervals of the 
distributions we calculate, assuming the tidal response of the s econdaries is domi nated by material strength. For 
Charon our theory is consistent with the numerical simulations of IStern et al.l ()2003l ). predicting an eccentricity about 
3 order of magnitudes smaller than observed. However, our analysis alleviates their nee d for numeri c al sim ulations as 
well as predicts the entire distribution of the eccentricity. The distributions measured bv I Stern et al.l (120031 ) are not all 
correct as their model includes only impact parameters out to twice the semi-major axis. In their simulations where 
q — 3.5 and 4.0 this excludes the impacts that are most relevant over an eccentricity damping timescale. Our results 
show that for q = 3.5 the interactions that dominate Charon's eccentricity are Pluto-sized perturbers interacting at 
about 200 times the semi-major axis! 

Even without eccentricity dissipation through tides, perturbations from other Kuiper belt objects are too weak to 
excite eccentricities of order 1 or inclination changes of order a radian for binaries that have orbital periods of a few 
days or weeks. It is not likely that the orbital planes of the close binaries have been affected significantly by other 
Kuiper belt objects given our current understanding of the history of the Kuiper belt. It falls on theories of binary 
formation to explain the distribution of orbital inclinations relative to the ecliptic for close binaries. Since edt) grows 
faster for binaries with large orbital periods, it is plausible that the smaller wide binaries (I998WW33 for example) 
have been brought to large eccentricities and inclinations by interacting with the rest of the Kuiper belt. 

When many binaries share the same perturbing swarm, such as in the Kuiper belt, we can use the eccentricities of 
all the binaries to probe the properties of the entire system. For example, if the mass spectrum is steeper than q = 4, 
the distribution of eccentricity is directly related to the slope and normalization of the mass spectrum. Conversely, 
the observed eccentricity can be used to place limits on the damping timescale of a binary and therefore the rigidity of 
those bodies. The small sample of Kuiper belt binaries with well measured eccentricities limits the current effectiveness 
of su ch a calculation. H owever, the Pan-STARRS project plans to detect around 20000 more members of the Kuiper 
belt (jl aiser et al]l2002l ): from these the number of orbit-determined Kuiper belt binaries will surely increase. 

The distribution we describe with equation [13] is a special case of a Levy distribution (Sato 199!^). This class 
of functions arise in the generalization of the central limit theorem to variables distributed with an infinite second 
moment. Alternatively, these functions can be characterized by the properties of the Levy flight they describe. For 
the eccentricity of the binaries discussed in this work, the frequency of a step is inversely proportional to a power of 
its size that depends on the mass spectrum of perturbers. It follows that the largest single step dominates the growth 
from accumulated smaller steps, causing, in the absence of damping, the typical eccentricity to grow faster than in 
a normal diffusive random walk. The slope of the distribution of excitations dictates the shape of the distribution. 
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This explains the coincidence of the distribution we derive in this work being exactly that of the distribution of 
eccentricity of protoplanets in a shear-dominated planetesimal disk, where the probability of changing; the ec centricity 
of a protoplanet is inversely proportional to the size of that change (iCollins fc Saril[200a[Coriins et al.l[200l . 

The authors thank Dmitri Uzdensky and Scott Tremaine for valuable discussions. R.S. is a Packard Fellow and an 
Alfred P. Sloan Fellow. This research was partially supported by the ERC. 

APPENDIX 

To calculate the excitation rates presented in sections |3] and [H it is necessary to integrate over all possible configu- 
rations of angles b and Vp relative to rf, and vt,. In this appendix we clarify the relation between the coefficients and 
equations [3] through [H 

We choose spherical polar coordinates for b and Vp to integrate equation [5) This requires a polar and azimuthal 
angle for b, 0}, and and a polar and azimuthal angle for Vp, 9y and 0^. By defining 9^ relative to b, the requirement 
that b and Vp be perpendicular fixes 9^ = 7r/2. 

The magnitude of the perturbation only depends on the vectors b and Vp relative to fh and Wb, so we use these 
vectors and their cross product, n to describe the components of &: h — hrVb + hyVi, + 6„n. The components are related 
to 9f, and 4>b in the typical way: hr — cos sin ^j, , h y — sin sin^b, and bn — cos We define the components of Vp 

relative to the same unit vectors. The angle (f>y describes the direction of Vp in the plane given by 6; the components 
of Vp follow from a rotation of this plane to align with n. We find the relations: 

Vr — bn COS (j)y — by{br sin — by cos (j)y) / {1 + bn) , 

Vv — bn sin + br{br sin (py — by cos (py) / {I + bn) , (1) 
Vn — ~br COS 4>y — by sin0^. 

The coefficient from equations ITT] and flTl (Ce), is defined to be the integral of |Ae|/(87r^(mp/mb)(wb/wp)(rb/6)^) as 
given by equation [T] 

(•^e) = / / / [{"^^rby + 2vrVy)^ + {\-vl- 2br)^Y'^ siii 9bd9bd(j)bd(f>y = 1.89 (2) 
47r Jq Jo Jo 

We similarly define (Ci) from equation [H 



r>27r /'27r ctt 

4^2 



(^«> = A / / / I'^^rbn + VrVnl sin 9bd9bd(f>bd(f,y = 0.75. (3) 

^ /o Jo Jo 

To calculate the coefficients used in the collisional excitation rate, equation[2Tl we use the |Ae| discussed in section 
[231 

{DJ-')^^!^ ^ j^{Avl+vl)^^-^^'^sm9bd9bd(pbd(py (4) 
For 7 = 2, the integral has a closed form solution of {D^) — E{—3), the complete Elliptic integral. For the inclination, 



(^r') = T^ / / \{vz)^"'\sin9bd9bd^bdc^y = - (5) 
47r^ Jo Jo Jo 7 

The coefficients for the excitation rates in the regime of 2 < 7 < 3 are more complicated as the dependence on b/rb 
cannot be factored out of the coefficient. In addition to integrating over all angles, we must integrate over impact 
parameter. For any 7, equation 1191 is written: 

pie) = z+i -7y — ^2-7^^2 , 6 

e^+^ \mb vq J zn 

where V.y-2 is discussed in section for a Gaussian distribution of perturber velocities, Vy-2 = 2r((l + 7)/2). 
The term {A2~^) again contains the angular information. Excitations for 2 < 7 < 3 are most important at 6 ^ r;, 
so we can not assume that b2 ~ b. We introduce explicit notation for the the components of the unit vector 
&2 — b2,rH + b2,yVb + b2.nn. Then the angular average coefficient is: 



-1 /'27r /'27T flT fOC 

'^'"'^^l ILL 



16 1 ^-.^V + 4^^'^'- 



,-, (7-l)/2 

xi sin 9bdxid9bd(f>bd(j)y, (7) 

X2 Xi / \ X2 Xl J 

with Xl = bjvb and X2 ~ b2/rb. The magnitude and components of b2 are related to b and Vp as described in section 

[2| b2 = b — Tb -I- Up (rf, • Vp). For 7 — 25/12 as discussed in l4.3[ (Ag'^^^^) w 15. For other 7 between 2 and 3, this factor 
is of the same order, 10-15. 
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